rm(list=ls(all=TRUE))
library(foreign)
library(sandwich)
library(lmtest)
library(metafor)

setwd("~/Dropbox/Apps/ShareLaTeX/Minimal Effects/replication/data")

master.sheet <- read.csv("master_sheet_input.csv", stringsAsFactors = FALSE)
master.sheet$label <- paste0(master.sheet$ExperimentName, ' - ', master.sheet$Seat)

##################
#HELPER FUNCTIONS#
##################

# Function to compute clustered standard errors, from Mahmood Arai.
cl <- function(fm, cluster){
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm), 2, function(x) tapply(x, cluster, sum))
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL)
}

analyze.data <- function(data, dv.t1, xvars, t0.range) {
  data <- read.dta(data)
  data <- subset(data, data[,"respondent_t1"] == 1 & data[,"t0_pid"] %in% t0.range)
  data.dv <- data[,dv.t1]
  # Rescale to mean 0 sd 1 in placebo group; treatment effects can then be interpreted
  # as the effect in standard deviations the treatment would have among an untreated population.
  data.dv <- (data.dv - mean(data.dv[!data$treat_ind], na.rm=TRUE)) /
    sd(data.dv[!data$treat_ind], na.rm=TRUE)
  xvars <- gsub("t0_pid ", "", xvars)
  t0.covariates.name <- unlist(strsplit(xvars, " "))
  x <- data[,t0.covariates.name]
  x <- as.matrix(x, dimnames = list(NULL, names(x)))
  
  lm.obj.covars <- lm(data.dv ~ data$treat_ind + x)
  result.covars <- cl(lm.obj.covars, data$hh_id)["data$treat_ind",]
  return(result.covars)
}
  
foo <- analyze.data(data = master.sheet[3,]$Data, 
                           dv.t1 = master.sheet[3,]$dv_t1, 
                           xvars = master.sheet[3,]$xvars,
                           t0.range = c(-3, -2, -1))

plot.data <- function(t0.range, title) {
  results <- data.frame(matrix(NA, nrow = 11, ncol = 5))
  names(results) <- c("Label", "Estimate", "SE", "t", "p")
  results$Label <- master.sheet[1:11, "label"]
  for(i in 1:11) {
    results[i, 2:5] <- analyze.data(data = master.sheet[i,]$Data, 
                                               dv.t1 = master.sheet[i,]$dv_t1, 
                                               xvars = master.sheet[i,]$xvars,
                                               t0.range = t0.range)
  }
  results$Estimate <- results$Estimate * 100
  results$SE <- results$SE * 100
  results$V <- results$SE^2
  meta <- rma.mv(yi = Estimate, V = V, data = results)
  pdf(paste0('../figures/',title,'.pdf'), width = 10, height = 4)
  forest(meta, cex = .8, slab = results$Label, xlab = paste0("Estimated Treatment Effect (CACE) in Percentage Points and 95% Confidence Interval Among ", title))
  dev.off()
}

plot.data(t0.range = c(-3, -2), "Republicans")
plot.data(t0.range = c(3, 2), "Democrats")
plot.data(t0.range = c(0, 1, -1), "Independents and Leaners")
plot.data(t0.range = c(-3, -2, -1), "Republicans and R Leaners")
plot.data(t0.range = c(3, 2, 1), "Democrats and D Leaners")
plot.data(t0.range = c(0), "Independents")
